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ABSTRACT 


In  this  thesis  we  present  an  algorithm  for  the  precise  determination  of  the  Mie 
extinction  efficiency  parameter.  The  mathematical  representation  of  the  Mie  parameters  is 
in  the  form  of  an  infinite  series,  and  any  technique  that  could  be  foimd  to  accelerate  the 
convergence  of  the  Mie  series  would  have  great  commercial  and  military  application. 
Results  are  presented  that  show  the  comparison  of  the  rate  of  convergence  obtained  by 
directly  summing  the  individual  terms  of  the  extinction  efficiency  parameter  and  the  rate 
obtained  using  an  existing  series  acceleration  technique.  It  was  found  that  the  acceleration 
method  we  employed,  known  as  the  Levin  method  of  series  transformation,  proved 
unsuccessful  in  accelerating  the  convergence  ofthe  Mie  series.  However,  other 
acceleration  techniques  exist  and  should  be  explored. 
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I.  INTRODUCTION 


Mie  scattering  is  the  scattering  of  electromagnetic  radiation  by  primarily  spherical 
particles  whose  diameters  are  comparable  to  the  wavelength  of  the  incident  radiation 
[Ref  1 :  p.  303].  In  nature,  we  see  this  manifested  as  the  white  appearance  of  clouds.  Due 
to  the  distribution  of  particle  sizes  in  the  clouds,  the  cumulative  effect  of  the  Mie  scattered 
Ught  waves  on  the  water  droplets  is  that  the  whole  spectrum  of  scattered  radiation 
combines  to  make  the  cloud  appear  white.  In  other  scenarios  which  involve  the  direct 
\ise  of  modem  technology,  Mie  scattering  has  many  applications  where  it  is  of  vital 
importance  to  have  a  highly  efficient  and  reliable  algorithm  for  signal  processing  and  for 
analyzing  the  properties  of  Mie  scattered  waves.  A  few  examples  include  the  reflection  of 
a  radar  signal  off  of  a  cloud  of  dispersed  particles  and  the  absorption  and  scattering  of 
radiation  from  soot  particles  in  the  atmosphere  [Ref  2];  Mie  scattering  as  a  technique  for 
the  sizing  of  air  bubbles  [Ref  3];  and  for  remote  sensing  applications.  Further  applications 
can  be  found  in  fields  as  diverse  as  astrophysies  [Refs.  4  -  5],  physical  chemistry  [Refs.  6- 
7],  and  a  unique  style  of  painting  that  employs  the  Mie  scattering  phenomenon  [Ref  8]. 

In  order  to  determine  the  amplitude,  extinction  efficiency,  or  scattering  efficiency 
of  an  electromagnetic  wave  that  has  undergone  Mie  scattering,  it  is  necessary  to  compute 
the  mathematical  representation  of  these  quantities,  or  Mie  parameters,  each  of  which  is 
analytically  represented  by  an  infinite  sum.  For  the  purpose  of  this  thesis  research,  our 
attention  has  been  focused  on  calculating  the  Mie  parameter  known  as  the  extinction 
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eflSciency  factor,  or  Qext,  which  describes  the  total  effect  of  scattering  and  absorption  in 
removing  radiation  from  the  incident  beam. 

As  will  be  shown  later  in  Chapter  IV,  the  number  of  terms  required  for  the  infinite 
series  to  converge  numerically  to  a  given  accuracy  is  directly  proportional  to  the  ratio 
between  the  wavelength  of  the  incident  radiation  and  the  circumference  of  the  scattering 
particle.  The  thrust  of  this  thesis  research  has  been  to  determine  if  this  infinite  series  can 
be  caused  to  converge  in  dramatically  fewer  iterations  through  the  utilization  of  a  series 
acceleration  technique.  The  basic  idea  is  to  transform  the  Mie  series  into  another  series 
that  converges  faster.  While  numerous  series  acceleration  methods  are  known,  one  of  the 
most  powerful  is  known  as  the  Levin  method  [Ref.  9:  p.  35].  The  ultimate  goal  of  this 
thesis  research  has  been  to  determine  if  the  Levin  method  can  significantly  accelerate  the 
convergence  of  the  Mie  series. 

As  an  example  of  the  power  of  the  Levin  method  (detailed  in  chapter  3),  consider 
the  alternating  harmonic  series  that  represents  the  natural  logarithm  of  2,  namely. 


n=\ 


(1) 


The  convergence  of  this  series  is  extremely  slow.  By  summing  the  individual  terms  of 
Eq.  (1),  over  a  million  and  a  half  terms  must  be  taken  to  achieve  an  accuracy  of  sfac  digits: 
1,565,239  terms  to  be  exact.  The  C  code  used  to  determine  this  value  and  the  program 
output  are  included  as  Appendices  A  and  B,  respectively.  By  emplojdng  the  Levin 
method,  however,  the  same  degree  of  numerical  accuracy  can  be  reached  after  only  sfac 
iterations!  The  ANSI  C  code  for  the  Levin  method  operating  on  the  alternating  harmonic 
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series  is  ineluded  as  Appendix  C.  After  the  fifteenth  term  of  the  Levin  method  is  reaehed, 
the  aceepted  value  of  aecuracy  out  to  the  sixteenth  digit  is  obtained  [Ref.  10:  p.  1 13]. 

This  is  an  astounding  reduetion  of  computer  processing  time  and  is  an  amaying 
demonstration  of  the  power  of  the  Levin  method  of  series  acceleration.  (We  note  that,  to 
obtain  eight  digit  accuracy  by  straightforward  summation  of  the  series  would  require  over 
one  hundred  million  terms;  nine  digit  accuracy  was  effectively  unattainable  using  the 
UNIX  system  resources  available  for  this  research.)  Table  1  shows  a  comparison  of  the 
partial  sum  obtained  fi-om  Eq.  (1)  and  those  calculated  with  the  Levin  method.  As  a 
further  comparison,  the  plot  shown  in  Fig.  1  was  produced  to  show  the  relative  speed  (i.e., 
number  of  iterations  required)  with  which  the  Levin  method  and  the  brute  force  method  of 
summing  the  individual  terms  of  the  alternating  harmonic  series  converge  to  99 


ln2  =  0.6931471805599453 


Iterations 

Partial  Sum 

Levin  Result 

1 

1.000000 

1.0000000000000000 

2 

0.500000 

0.6666666666666666 

3 

0.833333 

0.6944444444444444 

4 

0.583333 

0.6931372549019608 

5 

0.783333 

0.6931439393939393 

6 

0.616667 

0.6931474019283136 

7 

0.759524 

0.6931471777900349 

8 

0.634524 

0.6931471800150043 

9 

0.745635 

0.6931471806012293 

10 

0.645635 

0.6931471805592415 

11 

0.736544 

0.6931471805598542 

12 

0.653211 

0.6931471805599532 

13 

0.730134 

0.6931471805599452 

14 

0.658705 

0.6931471805599454 

15 

0.725372 

0.6931471805599453 

Table  1 .  Levin  Method  on  Logarithm  of  2. 
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Partial  Sum  /  In2 


Figure  1.  Levin  Method  vs.  Infinite  Series  for  ln(2) 
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percent  of  the  correct  value  of  ln(2).  Appendices  D  through  F  show  the  pertinent  ANSI  C 
and  Matlab  code.  Note  that  the  Levin  method  required  only  three  terms  to  reach  99 
percent  of  the  value  of  ln(2),  while  summing  the  individual  terms  of  the  series  required 
sixty-five  terms.  As  a  final  comparison,  the  time  required  to  reach  six  digits  of  accuracy 
using  the  brute  force  method  was  determined  (Appendices  A  and  B),  and  the  code  used  in 
Appendix  C  was  altered  to  determine  the  number  of  times  the  Levin  method  could  reach 
the  same  degree  of  accuracy  in  the  same  time  span  (1 1.7  seconds).  The  value  obtained 
was  approximately  20,000,  indicating  that  in  this  instance  the  Levin  method  is  about 
20,000  times  faster  than  summing  the  individual  terms  of  the  infinite  series. 

In  light  of  the  preceding  example,  the  question  arises  as  to  whether  the  Levin 
method  can  accelerate  the  Mie  series.  A  computer  program  was  written  to  calculate  the 
individual  terms  of  the  infinite  series  representation  of  the  Mie  parameter  Qext  and  to 
determine  the  number  of  these  terms  required  for  the  series  to  converge  to  an  accuracy  of 
one  part  in  10®.  Next,  the  Levin  method  was  applied  to  this  infinite  series,  and  the  number 
of  iterations  of  the  Levin  method  required  to  achieve  a  similar  accuracy  was  likewise 
determined.  A  comparison  of  these  two  results  showed  that  the  Levin  method  did  not 
succeed  in  accelerating  the  convergence  of  the  Mie  series.  In  fact,  the  number  of 
iterations  required  for  the  Levin  method  to  converge  was  comparable  to  that  obtained  by 
simply  summing  the  infinite  series  representation  of  the  Mie  parameter.  Considering  the 
algorithm  which  the  Levin  method  utilizes,  this  method  when  apphed  to  the  Mie  series  will 
actually  result  in  increased  processor  time.  Although  the  Levin  method  proved 
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unsuccessful  in  accelerating  the  convergence  of  the  Mie  series,  further  research  in  this  area 
is  still  needed  as  there  are  other  series  acceleration  methods  worth  exploring 
[Ref  11:  p.  56]. 

In  the  next  section,  we  discuss  the  relevant  Mie  parameters  and  the  infinite  series 
for  computing  the  extinction  coefficient.  In  Chapter  III  we  discuss  the  Levin  convergence 
acceleration  method.  In  Chapter  IV  we  discuss  our  results  for  applying  the  Levin  method 
to  the  Mie  series. 
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II.  THE  MIE  SCATTERING  PARAMETERS 


A.  BACKGROUND 

There  are  four  basic  Mie  scattering  parameters.  The  first  two  pertain  to  the 
scattering  amplitude  and  describe  the  complex  amplitudes  of  the  perpendicular  and  parallel 
components  (with  respect  to  the  scattering  plane)  of  the  electric  field  vector.  The  third 
Mie  parameter  is  the  scattering  efficiency  factor,  while  the  fourth  parameter  is  the 
extinction  efficiency  factor.  All  fom  of  these  expressions  are  given  in  dimensionless  terms, 
and  the  evaluation  of  these  four  basic  functions  rehes  on  the  proper  evaluation  of  four 
coefficients:  the  Mie  coefficients  a„  and  b„,  and  the  angular  coefficients  n„  and  r„. 

These  four  expressions  can  be  described  using  three  basic  parameters:  a 
dimensionless  size  x  of  the  scattering  sphere,  an  index  of  refi-action  m  (usually  complex), 

and  a  scattering  angle  ^(relative  to  the  forward  direction).  The  dimensionless  size 
parameter  x  is  given  by 

^  =  ^  circumference  of  sphere 

\  ■^ )  wavelength  ’ 

where  k  =  lnl?L  is  the  wave  number,  r  is  the  radius  of  the  sphere,  and  A  is  the 
wavelength  of  the  incident  radiation.  The  complex  index  of  refraction  m  is  given  by 

m=v-Ki  ,  (3) 

where  v  is  the  real  part  of  the  index  of  refi-action  and  k  is  the  imaginary  part.  A  complex 
refi-active  index  indicates  an  absorbing  sphere,  such  as  soot.  A  purely  real  index  of 
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refraction  is  indicative  of  a  non-absorbing  sphere;  for  the  purpose  of  this  thesis  research, 
only  cases  of  purely  real  refractive  index  were  considered. 

B.  SCATTERING  AMPLITUDE  PARAMETERS 

Let  the  electric  field  vector  amplitude  of  the  scattered  radiation  field  be  given  by 
Asc.  This  radiation  field  can  in  turn  be  expressed  in  terms  of  the  scalar  perpendicular  and 
parallel  components  Ai  and  Az.  We  can  define  dimensionless,  complex  amplitudes  Si  and 
52  by  multiplying  the  amplitudes  and  Az  by  the  free-space  propagation  constant  k, 
which  in  turn  can  be  represented  by  an  infinite  converging  series  [Ref  1 2 :  p.  1 3] : 


The  Mie  coefficients  and  b„  are  fimctions  of  the  index  of  refraction  m  and  the 
size  X,  while  the  angular  coefficients  ;r„  and  t„  are  fimctions  of  ju=cos9  only.  The  latter 
coefficients  are  defined  in  terms  of  Legendre  polynomials  and  their  derivatives.  The 
coefficients  a„  and  b„  are  functions  of  spherical  Bessel  functions  of  the  first,  second,  and 
third  kind,  and  can  be  expressed  in  a  variety  of  ways. 

C.  SCATTERING  EFFICIENCY  AND  EXTINCTION  PARAMETERS 

It  can  be  shown  that  the  differential  scattering  cross  section  for  unit  incident  flux  is 
dcx{m,  x,0)  =  ^  Age  •  Alc(m,x,ff)d(o  ,  (5) 


where  dco  is  an  element  of  the  solid  angle. 
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Following  Diermendjian  [Ref.  12:  p.  13],  by  representing  the  unpolarized  incident 
radiation  as  the  sum  of  two  independent  and  linearly  polarized  components  of  equal  flux, 
the  scattering  cross  section  can  be  expressed  as 

=  J  d(T(m,x,0)  =  ^  J(  A4*  +  AA*)do)  ,  (6) 

£1  ^  n 

where  Q  =  4;r  is  the  solid  angle. 

The  scattering  eflSciency  factor  Qg(.(m,x)  is  obtained  by  normalizing  a by  the 
sphere’s  geometrical  cross  section  : 

^  ^  n=I 

The  total  extinction  cross  section  and  efficiency  fector,  which  includes  the 
contribution  due  to  absorption,  can  be  similarly  defined.  This  leads  to  the  following 
expression  for  extinction  efficiency  factor: 

2'exr=:4E(2"  +  l)M^»+^«)  •  (8) 

X 

As  noted  earlier  for  the  case  of  a  purely  real  index  of  refraction  m,  there  will  be  no 
absorption,  and  thus  the  extinction  efficiency  factor  Qext  will  be  identically  equal  to  the 
scattering  efficiency  factor  Qsca  . 

D.  THE  MIE  COEFFICIENTS 


The  simplest  and  most  elegant  expressions  for  the  Mie  coefficients  are  as  follows 
[Ref  13:  p.  195]: 


J„(kr) 

h^\ilcr)  ’ 


(9a) 
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(9b) 


ly  ^JAf^)-fVn(kr) 
krh^_^{kr)  -  nh^'\{kr) 

where  jn(kr)  is  a  spherical  Bessel  function  of  the  first  kind  with  order  n,  U‘^„(kr)  is  a  Bessel 
function  of  the  third  kind,  known  as  a  Hankel  function,  and  is  defined  as 

h^^\(kr)  =  j„(kr)  +  iy„(kr)  ,  (10) 

and  y„(kr)  is  a  spherical  Bessel  function  of  the  second  kind  with  order  n. 

Another  common  form  of  the  Mie  coefficients  is  that  adopted  by  van  de  Hulst 
[Ref  14:  p.  123]: 


(11a) 


(lib) 


where  z  =  mx ,  while  and  are  Ricatti-Bessel  functions  which  are  defined  as 
follows: 


Wn{x)  =  Xj„(,x)  , 

(12a) 

X„ix)  =  -xy„{x)  , 

(12b) 

=  .  (12c) 

Here,y„  and^„  are  spherical  Bessel  functions  of  the  first  and  second  kind,  respectively. 

For  computational  purposes,  it  is  best  to  express  the  Mie  coefficients  in  a  form 
conducive  to  separating  their  real  and  imaginary  components.  Through  the  use  of 
recursion  formulas  and  circular  functions  [Ref  12:  p.  16-19],  the  Mie  coefficients  may  be 
written  in  the  alternative  form  [Ref  15:  pp.  16-17]: 
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\  m  x] 

d  =  - 1 — 

n  (a 


^+^[CW-C,W 


(13a) 


.J 


n\ 


-  y/„-^{x) 

l'”A(^)+^|c(^)-Ci(^) 


(13b) 


These  forms  of  the  Mie  coefiBcients  were  programmed,  using  the  C  language,  to  determine 
the  Mie  coefiBcients  and  hence  the  extinction  efiBciency  factor  Qext . 

E.  CALCULATING  USING  ANSI  C 


The  author  wishes  to  thank  the  writers  of  "Numerical  Recipes  in  C"  [Ref.  16]  for 
enabling  the  circumvention  of  the  reinvention  of  the  wheel,  in  that  their  code  for 
computing  ordinary  and  spherical  Bessel  functions  and  their  derivatives  was  utilized  in  this 
thesis  research  for  computing  the  Mie  coefiBcients  and  calculating  the  extinction  efiBciency 
factor.  Of  noteworthy  interest  is  that  said  code  cites  Jerry  Lentz  of  the  Naval 
Postgraduate  School  as  the  author  of  an  inqjroved  technique  for  calculating  Bessel 
functions  through  the  use  of  continued  fractions. 

The  code  used  for  calculating  the  ordinary  Bessel  functions  is  incorporated  into  a 
structure  which  returns  the  values  of  the  Bessel  functions  of  the  fast  and  second  kind,  and 
their  fibrst  derivatives.  The  spherical  Bessel  ftmction  code  simply  generates  a  normalization 
factor,  makes  a  call  to  the  function  which  calculates  ordinary  Bessel  functions  (of  half-odd 
integer  order),  then  in  like  fashion  returns  values  of  the  spherical  Bessel  functions  of  the 
first  and  second  kind,  and  their  first  derivatives.  The  code  for  calculating  the  ordinary 
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Bessel  functions  also  makes  calls  to  two  additional  functions,  "beschb"  and  "chebev",  also 
obtained  from  "Numerical  Recipes  in  C"  [Ref.  16].  The  “beschb”  function  is  used  to 
evaluate  the  gamma  functions  present  in  the  Bessel  function  formula,  calling  the  “chebev” 
function  in  the  process  to  perform  Chebyshev  expansion. 

Once  the  interfacing  codes  for  calculating  ordinary  and  spherical  Bessel  functions 
were  brought  together,  a  program  was  written  to  assemble  the  components  of  the  Mie 
coefficients  as  expressed  in  Eqs.  (i3a)  and  (13b).  Again,  the  Mie  coefficients  are 
fimctions  of  the  dimensionless  size  parameter  x  and  the  index  of  refraction  m. 

Next,  a  loop  was  created  that  calculated  the  individual  terms  of  the  extinction 
efficiency  factor  (one  of  the  Mie  parameters),  each  term  of  which  relied  on  a  calculated 
value  of  the  Mie  coefficients.  As  is  the  case  of  all  of  the  Mie  parameters,  the  extinction 
efficiency  factor  is  represented  as  an  infinite  sum.  Thus,  the  loop  for  calculating  the 
efficiency  extinction  feictor  was  summed  xmtil  the  value  converged  to  one  part  in  10*. 

With  each  increment  of  the  loop  counter,  the  Bessel  functions  associated  with  the  Mie 
coefficients  were  of  subsequent  higher  order,  i.e.,  for  n=l,  J„ix)  would  be  a  Bessel 
function  of  the  first  order;  for  n=2,  J„{x)  would  be  a  Bessel  function  of  the  second  order, 
and  so  on. 

As  a  test  for  the  accuracy  of  the  code  and  the  correctness  of  the  programming,  a 
routine  was  created  that  calculated  the  extinction  efficiency  factor  for  size  parameter  x 
varying  from  0  to  9,  using  a  step  size  of  0.01  and  applying  6  different  values  for  the 
refractive  index.  A  plot  showing  all  six  curves  was  produced  from  the  data  obtained  from 
this  program  and  is  shown  in  Fig.  2.  This  plot  very  closely  resembles  that  which  was 
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produced  by  Van  de  Hulst  [Ref.  14:  p.  151],  with  the  exception  that  Van  de  Hulst's  plot 
was  intentionally  smoother  as  a  result  of  not  showing  fine  detail.  The  C  code  used  to 
produce  data  for  Fig.  2  is  included  as  Appendix  G,  with  the  output  data  being  sent  to  a  file 
called  “Qout”.  The  Matlab  code  used  to  produce  Figure  2  is  included  as  Appendix  H. 
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III.  THE  LEVIN  METHOD 

A.  SERIES  ACCELERATION  METHODS 

The  main  idea  behind  the  Levin  method,  as  in  any  series  acceleration  technique,  is 
to  transform  a  given  slowly  converging  series  into  another  series  that  converges  faster. 

We  note  here  the  distinction  between  a  sequence  and  a  series.  An  infinite  sequence  (5«)  is 
an  unending  progression  of  numbers  S„  which  may  be  real  or  complex.  A  sequence 
converges  if  a  number,  S,  exists  so  that,  corresponding  to  every  positive  number  s,  no 
matter  how  small,  a  number  can  be  found  such  that  \S„  —  S’]  <  £  for  n>n^.  In  this  case, 

the  sequence  (S„)  is  said  to  converge  to  the  limit  S  as  n  tends  to  infinity.  An  infinite  series, 
on  the  other  hand,  is  the  sum  of  an  infinite  sequence.  Let  , . . . ,  , . . .  be  an  infinite 

sequence  of  numbers,  real  or  complex.  Let  the  sum  «,  +  M2  ^  denoted  by  S„, 

which  is  called  the  partial  siun.  Then,  if  the  sequence  of  partial  sums  (S„)  converges  to 
a  limit,  S,  the  infinite  series  m,  +  M2  h —  is  said  to  be  convergent,  or  to  converge  to  the 
sum  S.  The  connection  with  the  Levin  method  is  to  transform  the  sequence  of  partial 
sums  into  another  sequence  that  converges  faster,  i.e.,  requires  fewer  terms  to  converge  to 
the  limit. 

Given  a  sequence  of  real  or  complex  numbers  (S„)  which  converges  to  S,  this 
sequence  can  be  transformed  into  another  sequence  (T„).  A  trivial  example  of  such  a 
sequence  transformation  is 

n  =  0,l,.,.  .  (14) 
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To  be  useful,  however,  the  following  properties  for  the  transformed  sequence  must  hold 
[Ref  11:  p.  1]: 

1 .  (T„)  must  converge, 

2.  (Tn)  must  converge  to  the  same  limit  as  (S„),  and 

3 .  (T„)  must  converge  to  S'  faster  than  (S„);  that  is,  lim (7^  -  S) /(S„  -  S)  =  0. 


If  property  3  holds,  the  transformation  T  is  said  to  accelerate  the  convergence  of  the 
sequence  (Sn),  or  that  the  sequence  {T„)  converges  faster  than  {S„).  An  example  of  a 
sequence  transformation  meeting  these  three  conditions  is  Aitken’s  process  [Ref  1 1 ; 
pp.  1-7]: 


T„=^^‘S„  + 
”  A'S„  " 


^  AC  ^ 
1  _  ^"+1 


V 


A"S, 


-’n+l  ’ 


(15) 


or. 


Tn=- 

S. 


Sn^n+2  ~ 


n+\ 


n+2 


•2S„,.+S„ 


n  =  0,l 


(16) 


where  A  is  the  difference  operator  defined  by  Av„  =  v„^.,  -  v„  and  A**'v„  =  A*v„^,  -  A*v„ . 


The  A^  in  the  denominator  of  Eq.  (15)  accounts  for  the  process’s  name. 

The  Levin  transformation,  as  with  most  other  transformation  algorithms,  is  a 
particular  case  of  what  is  known  as  the  E-transformation,  which  is  the  most  general 
sequence  transformation.  Transformations  belonging  to  the  E-transformation  class  include 
[Ref  1 1:  p.  56]:  Richardson  polynomial  extrapolation.  Shanks’  transformation,  the  Levin 
transformation,  and  many  other  historically  known  acceleration  methods.  We  note  that 
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series  acceleration  is  an  active  field  of  research  in  numerical  analysis,  with  the  Levin 
method  and  E-transformation  having  been  discovered  only  in  the  past  25  years. 

B.  THE  E-ALGORITHM 


Following  Brezinski  and  Zaglia  [Ref  1 1:  pp.  56-57],  the  E-transformation  is  based 
on  the  following  relation: 

S„-S-a^g,(n) - a4g^(«)  =  0  ,  (17) 

where  5n  is  an  element  of  the  sequence  to  be  transformed,  the  a’s  are  unknown  scalars,  the 
g,(n)’s  are  given  auxiliary  sequences  (which  may  depend  on  terms  of  the  sequence  S„ 
itself),  and  where  k  is  a  fixed  integer.  Rewriting  Eq.  (17)  in  terms  of  iS„, 

S„  =  S  +  a,g,(n)  +  ---  +  a„g^(n)  .  (18) 

The  basic  idea  behind  Eq.  (18)  is  to  attempt  to  fit  the  actual  behavior  of  S„,  as  a  function 
of  n,  so  that  it  may  be  extrapolated  smoothly  to  the  (unknown)  limit,  S.  This  fitting  is 
achieved  through  the  particular  choices  of  the  fimctions  gX”)-  By  incrementing  n  in  Eq. 
(18)  to  n  +  ^,  one  has  A:  +  1  equations  in  A  +  1  unknowns  (the  limit  S  and  the  scalars 
a, ,  •  •  • ,  ).  By  solving  these  equations,  we  obtain  a  sequence  of  estimates  for  the  limit. 


which  we  denote  by 


- 

gi(n)  ••• 

g,(n  +  A) 

gt(n)  "• 

8kif^  +  k) 

1 

1 

&(«)  ••• 

8i(ri  +  k) 

gt(«)  ••• 

gt(.n  +  k) 

(19) 
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It  is  assumed  that  the  determinant  in  the  denominator  of  Eq.  (19)  is  not  equal  to  zero. 
Finally,  note  that  Eq.  (18)  is  the  kernel  of  the  transformation  Eq.  (19);  that  is,  Eq.  (18)  is 
the  set  of  sequences  for  which  there  exists  the  sequence  S  such  that,  for  all  n,  =  S. 

The  E-algorithm  is  the  recursive  algorithm  that  allows  one  to  compute  the  numbers 
without  actually  computing  the  determinants  in  Eq.  (19).  Given  the  following  rules. 


£<">  =  S, , 


«o?=a(").  B  =  0,1,...  and  1  =  1,2,...  , 


then  for  A:  =  1,2,...  and  n  =  0,  1, . . . ,  the  main  rule  for  the  E-algorithm  is 


£<■>  =  £<!>  - 


^k-\  ^k-\ 


„(«+!) 

Sk-\,k 


.p(”) 

6i-U- 


(20a) 

(20b) 


(21) 


where  the  auxiliary  sequences  computed  by  the  following  auxiliary  rule: 


(n+l)  _  («) 

o-W  _  («) 

&k-\,i  _  fSn)  ok-\,k 

Sk-\,k  Sk-\,k 


i  =  k  +  \,k  +  2,... 


(22) 


C.  LEVIN’S  TRANSFORMS 


Levin’s  method  of  generating  non-linear  transformations  for  increasing  the  rate  of 
convergence  of  sequences  [Ref  17:  pp.  371-388]  was  first  introduced  in  1973.  For 
Levin’s  generalized  transformation,  the  auxiliary  sequences  denoted  by  the  g/(n)’s  in  the 
E-transformation  given  in  Eq.  (19)  take  on  the  particular  form: 


giip)  = 


n _ ^ 


(23) 


where  x„  andy„  are  themselves  auxiliary  sequences. 
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Levin’s  transforms  are  basically  generalizations  of  Aitken’s  process  and  of  the 
E-transformation  corresponding  to  the  first  column  of  the  E-algorithm  as  shown  in 
Eq.  (21)  [Ref.  1 1 :  p.  1 13],  denoted  as  £, .  The  kernel  of  Aitken’s  process  is 

S„-S  =  a-AS„  ,  (24) 

while  the  kernel  of  the  transformation  is  of  the  form 

S„-S  =  t„-g{n)  .  (25) 

By  defimtion  in  the  Levin  method,  g(ri)  is  taken  to  be  an  arbitrary  polynomial  of  degree 

—  l)  of  the  quantity  («)  .  Also,  tn  denotes  the  n**  element  of  the  series  (as  for  example 
the  term  of  the  Mie  series).  Thus,  the  sequences  of  Eq.  (25)  can  be  expressed  as 

S„-S  =  t„- a^+a2{n)  +*--  +  aj(«)  .  (26) 

FoDowing  Brezinski  and  Zaglia,  multiplying  both  sides  of  Eq.  (26)  by  (w)*''  and 


rearranging  terms  yields 


/  x(*-l)  s  ~S  ,  at-l)  ,  ./t~2 

(«)  -^  =  a,(«)'  ’+a2(n)  +-+a,  . 


By  applying  the  operator  A*  to  both  sides  of  Eq.  (27),  the  right  hand  side  becomes 
identicaUy  zero,  since  it  is  a  polynomial  of  degree  (A:  - 1)  in  n.  Therefore,  for  all  n. 


\  Li 


Since  A*  is  a  linear  operator,  the  above  expression  can  be  also  be  expressed  as 


r  L  t. 
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The  numbers  ,  denoted  now  as  the  Levin  estimate  ,  are  therefore  given  by 


(30) 


Using  the  well-known  formula  [Ref  1 1 :  p.  115] 


A"". =2;(-i)‘ ck*. . 


(31) 


where  is  the  binomial  coefficient  defined  by 


Nl 


k\{N-k)\  ’ 


Eq.  (30)  then  becomes 


(«)  -izL 


2](-l)*(n  +  A:)"‘'^iV!  lk\{N-k)\ 


VC  = 


*«+A: 


'N  ^  N 


k=\  ^n+k  I 


(32) 


(33) 


Equation  (33)  defines  the  sequence  of  Levin  estimates  for  the  limit  of  the  series. 

D,  ENCODING  THE  LEVIN  TRANSFORM 


In  the  instances  where  the  Levin  code  (Appendices  C  and  L)  was  utilized,  TV  was 
used  to  represent  the  partial  sum  fi'om  Eq.  (33)  above.  The  sum  TV  is  related  to  the 
sequence  4  by  the  following  expression: 

T«  =  14  .  (34) 

i=1 

In  other  words,  4  represents  the  individual  term  for  the  loop  which  calculates  the 
extinction  efficiency  fector  ,  and  TV  is  the  value  for  which  converges  to  one  part 
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in  10*.  Equation  (33)  was  split  into  two  individually  calculated  terms,  and 

representing  the  numerator  and  denominator,  respectively.  Consequently,  the  following 
equations  result: 

=  — ,  and  (35) 

*=i  h  / 

e[A']  =  i;(-l)‘W''''^Af!A!(JV-t)!  .  (36) 

*=i  h  / 

Thus,  Eq.  (30)  becomes 

Equation  (34)  was  encoded  in  order  to  calculate  the  Levin  estimate  of  ,  and  the 
number  of  terms  needed  for  the  Levin  transform  to  cause  the  infinite  series  representation 
of  to  converge  was  subsequently  determined.  This  number  was  compared  with  that 
which  was  obtained  by  direct  summing  of  the  individual  terms  of  the  infinite  series 
representation  of  to  convergence  in  order  to  determine  if  the  Levin  method 
accelerated  the  convergence. 
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Figure  2.  Extinction  Efficiency  Factor  vs.  Size  Parameter 
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IV.  RESULTS 


After  completing  the  Mie  code  and  producing  values  for  the  extinction  efficiency 
factor  that  agreed  with  the  tabulated  results  of  Van  de  Hulst,  the  result  obtained  from 
applying  the  Levin  method  of  series  acceleration  to  the  iirhnite  series  representation  of  a 
Mie  scattering  parameter  was  somewhat  anti-climactic,  albeit  not  entirely  unanticipated. 
Simply  put,  the  Levin  method  did  not  accelerate  the  convergence  of  the  extinction 
efficiency  factor.  In  point  of  fact,  the  number  of  iterations  required  for  the  Levin  method 
to  converge  to  a  given  accuracy  was  approximately  the  same,  less  one  or  two  terms,  as  the 
number  of  terms  required  by  direct  summation.  Size  parameters  ranging  from  0  to  100 
were  used  in  making  this  determination.  Since  the  Levin  method  requires  calls  to 
additional  functions  that  calculate  the  log  of  a  gamma  function,  the  log  of  a  factorial,  and  a 
binomial  coefficient  at  each  iteration  of  the  summation  loop  -  the  codes  for  which  were 
obtained  from  “Numerical  Recipes  in  C”  [Ref  16]  -  this  method  is  actually  slower  than 
direct  summation.  It  should  also  be  noted  that,  for  size  parameters  x  >  108 ,  the  Levin 
method  fails  entirely  because  the  individual  terms  U  become  vanishmgly  small,  producing 
undefined  output.  The  code  used  to  determine  via  summing  the  infinite  series  and  the 
code  utilizing  the  Levin  estimate  are  included  as  Appendices  I  and  J,  respectively. 

The  C  code  used  to  produce  the  data  used  for  the  plot  shown  in  Figure  2  was 
slightly  altered  to  determme  the  number  of  iterations  required  for  the  Mie  series  to 
converge  as  a  function  of  size  parameter.  The  result  is  shown  in  Figure  3.  Since  both 
methods  of  convergence  required  the  same  number  of  iterations,  the  code  used  for 
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summing  the  individual  terms  of  the  Mie  series  (Appendix  G)  Avas  chosen  to  generate  the 
data  used  for  the  plot  shown  in  Figure  3.  A  linear  relationship  can  easily  be  seen  for  all 
indices  of  refraction,  with  the  number  of  terms  in  the  series  proportional  to  the  size 
parameter. 

While  a  tremendous  amount  of  effort  was  put  forth  to  come  to  the  conclusion  that 
the  Levin  method  is  unsuccessful  in  accelerating  the  convergence  of  the  Mie  series,  the 
work  was  not  all  for  naught.  An  important  question  was  answered,  leading  the  way  for 
further  research  in  this  area. 


Figure  3.  Rate  Convergence  vs.  Size  Parameter 
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V.  CONCLUSION 


In  this  thesis  research  we  assembled  the  code  to  calculate  the  Mie  coefficients, 
integral  components  of  all  four  of  the  Mie  scattering  parameters.  Next,  we  calculated  the 
Mie  extinction  efficiency  factor  by  summing  the  terms  of  the  infinite  series  representation 
to  convergence.  Several  different  indices  of  refraction  were  used,  and  a  plot  of  the 
extinction  efficiency  factor  versus  size  parameter  was  produced.  These  results  were  formd 
to  agree  closely  with  that  of  Van  de  Hulst’s  authoritative  work,  thereby  ensuring  the 
accuracy  of  the  code.  Finally,  an  algorithm  for  the  Levin  method  of  series  transformation 
was  incorporated  into  the  existing  code.  The  rate  of  convergence  was  subsequently 
determined  and  compared  with  the  rate  of  convergence  achieved  by  summing  the 
individual  terms  of  the  Mie  series.  Based  on  the  results,  we  arrived  at  several  important 
conclusions. 

First,  the  number  of  terms  required  to  converge  was  found  to  be  linearly  related  to 
the  size  parameter.  This  included  size  parameters  both  within  and  outside  the  Mie  regime. 
The  linear  relationship  continued  until  a  limitation  was  reached,  which  leads  us  to  our 
second  conclusion. 

The  magnitude  of  the  size  parameter  was  found  to  impose  a  limit  on  the  extent  to 
which  the  Levin  method  was  able  to  cause  the  Mie  series  to  converge.  While  a  size 
parameter  of  108  is  well  beyond  the  Mie  scattering  region  and  within  the  optical  region,  it 
stiU  points  to  a  computational  limit  of  the  Levin  method. 
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Most  importantly,  we  have  coneluded  that  the  Levin  method  of  series 
transformation  did  not  accelerate  the  convergence  of  the  Mie  series.  Several  simulations 
were  run  using  size  parameters  ranging  from  less  than  one  up  to  one  himdred,  and  indices 
of  refraction  ranging  from  one  to  two.  All  results  showed  a  rate  of  convergence,  insofar 
as  the  number  of  iterations  required,  equal  to  the  rate  obtained  by  summing  the  individual 
terms  of  the  Mie  series.  This  is  an  important  result.  There  are  other  series  transformation 
methods  in  existence  that  may  accelerate  the  convergence  of  the  Mie  series,  and  these 
should  be  explored  in  their  own  right. 
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APPENDIX  A:  ALTERNATING  HARMONIC  SERIES  CODE  IN  ANSI  C 


/*  Thesis  Program  to  calculate  number  of  terms  for  alternating  harmonic  series  to  reach  progressive  */ 
/*  digits  of  accuracy  */ 

/*  LT  Brian  Johnson  */ 

/*  Compiler:  Borland  C++  Ver.  5.0  */ 

/*  File  Name:  ln2x  */ 

#  include  <stdio.h> 

#  include  <math.h> 

#  include  <time.h> 

int  n; 

int  places  =  1; 

double  AltHarmSum  =  0.0; 
double  RoundedNew  =  0.0; 
double  RoundedOld  =  0.0; 

float  accuracy  =  10.0; 

main() 

{ 

n  =  0; 

do  { 
n+=l; 

AltHarmSum+=(double)(  pow(-l,n+l)/n); 

RoundedNew  =  floor(AltHarmSum*accuracy+.5)/accuracy; 

if  (  RoundedNew  ==  RoundedOld ) 

{ 

printfl["\nlt  took  %d  terms  to  reach  %d  digits  of  accuracy:”,  n,  places); 

printf("\n  Sum  =  %f ",  AltHarmSum); 

printfl;”\n  Time  elapsed:  %3.2f  sec.\n",  clock()/le6  ); 

places+=l; 

accuracy*=10.0; 

} 

RoundedOld  =  RoundedNew; 

}  while  (places  <=  15); 

} 
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APPENDIX  B.  OUTPUT  OF  ALTERNATING  HARMONIC  SERIES  CODE 

It  took  12  terms  to  reach  1  digits  of  accuracy: 

Sum  =  0.653211 
Time  elapsed:  0.00  sec. 

It  took  271  terms  to  reach  2  digits  of  accuracy: 

Sum  =  0.694989 
Time  elapsed:  0.00  sec. 

It  took  1417  terms  to  reach  3  digits  of  accuracy: 

Sum  =  0.693500 
Time  elapsed:  0.01  sec. 

It  took  177341  terms  to  reach  4  digits  of  accuracy: 

Sum  =  0.693150 
Time  elapsed:  1.33  sec. 

It  took  229300  terms  to  reach  5  digits  of  accuracy: 

Sum  =  0.693145 
Time  elapsed:  1.71  sec. 

It  took  1565239  terms  to  reach  6  digits  of  accuracy: 

Sum  =  0.693147 
Time  elapsed:  1 1 .70  sec. 
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APPENDIX  C:  LEVIN  ALGORITHM  ON  ALTERNATING  HARMONIC  SERIES 


/*  Thesis  Program  to  apply  Levin  method  to  alternating  harmonic  series  and  compute  number  of  terms  */ 
/*  for  series  to  converge  */ 

/*  LT  Brian  Johnson  */ 

/*  Compiler;  Borland  C++  Ver.  5.0  */ 

/*  File  Name:  ln21ev.c  */ 

#include  <stdio.h> 

#include  <stdlib.h> 

#include  <stddef  h> 

#include  <math.h> 

#define  MAX  15 

int  k; 

double  T[MAX+1],  t[MAX+l]; 

double  a; 

main() 

{ 

T[0]  =  t[0]  =  0.0; 
k=  1; 

printf("\n\ni  T[i]  t[i]  P[i]/Q[i]\n"); 

while(k<=MAX) 

{ 

tpc]  =  (doubleX  pow(-l,k+l)/k); 

T[k]  =  T[k-l]  +  t[kl; 
k  =  k  +  1; 

} 

a  =  levin(); 

} 

double  levinQ 

{ 

double  P[MAX+1],Q[MAX+1]; 
int  1; 

1=1; 

P[0]  =  Q[0]  =  0.0; 
while(l<=MAX) 

{ 

P[l]=P[l-l]+pow(-l,l)*pow(l,MAX-l)*T[l]/t[l]*bico(MAX,l); 

Q[1]=Q[1- 1  ]+pow(- 1  ,l)*pow(kMAX- 1  )/t[l]*bico(MAX,l); 

1=1+1; 

printfl:"\n%d  %f  %f  %f",i,T[i],t[i],P[i]/Q[i]); 

} 

return  (P[MAX]/Q[MAX]); 

} 
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double  bico(int  n,  int  k) 

{  double  factln(int  n); 

return  £loor(0.5+exp(factln(n)-factln(k)-factln(n-k))); 

} 

double  factln(int  n) 

{  double  gammlnCfloat  xx); 

static  float  a[101]; 
if  (n  <=  1)  return  0.0; 

if  (n  <=  100)  return  a[n]  ?  a[n] :  (a[n]=gammln(n+1.0)); 
else  return  ganunln(n+1.0); 

} 

double  ganunln(float  xx) 

{  double  x,y,tmp,ser, 

static  double  cofI6]={76.18009172947146,-86.50532032941677, 
24.0 140982408309 1,-1 .23 1739572450155, 
0.1208650973866179e-2,-0.5395239384953e-5}; 

intj; 

y=x=xx; 

tmp=x+5.5; 

tmp  —  (x+0.5)*log(tmp); 

ser=1.000000000190015; 

for  (j=0;j<=5y-H-)  ser  +=  coflj]/++y; 

return  -tmpt-log(2.50662827463 10005*ser/x); 

} 
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APPENDIX  D:  ANSI  C  CODE  TO  PRODUCE  DATA  FOR  ALTERNATING 

HARMONIC  SERIES  PLOT 


/*  Thesis  Program  to  produce  data  for  plot  of  convergence  of  alternating  harmonic  series  to  ln(2)  */ 
/*  LT  Brian  Johnson  */ 

/*  Compiler:  Borland  C++  Ver.  5.0  */ 

/*  File  Name:  ln2plot.c  */ 

#  include  <stdio.h> 

#  include  <math.h> 

#  include  <time.h> 

int  n; 

int  places  =  1; 

double  AltHarmSum  =  0.0; 
double  RoimdedNew  f  0.0; 
double  RoundedOld  =  0.0; 

float  accuracy  =  10.0; 

main() 

{ 

for  (n  =  1;  n  <=  51;  n++)  { 

AltHarmSum+=(double)(  pow(-l,n+l)/n); 

RoundedNew  =  floor(AltHarmSum*accuracy+.5)/accuracy; 

printf(”\n%d  %1.3f\  n,  log(2)  -  fabs(log(2)  -  AltHarmSum)); 

} 

} 
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APPENDIX  E.  ANSI  C  CODE  TO  PRODUCE  DATA  FOR  PLOT  OF  LEVIN 
METHOD  ON  ALTERNATING  HARMONIC  SERIES 

/*  Thesis  Program  to  produce  data  for  plot  of  Levin  method  on  the  convergence  of  alternating  */ 

/*  harmonic  series  to  ln(2)  */ 

/*  LT  Brian  Johnson  */ 

/*  Compiler:  Borland  C++  Ver.  5.0  */ 

/*  File  Name:  ln21evplot.c  */ 

#include  <stdio.h> 

#include  <stdlib.h> 

#include  <stddef  h> 

#include  <math.h> 

#include  "nrutil.h" 

#include  "nrutil.c" 

#include  <time.h> 

#defineMAX16 

double  levinQ; 
double  bico(int,  int); 
double  factln(int); 
double  gammln(float); 
double  T[MAX+1],  t[MAX+l]; 
double  a; 

int  k,  i; 

mainO 

{ 

T[0]  =  t[0]  =  0.0; 
k=l; 

while(k<=MAX) 

{ 

t[k]  =  (double)(  pow(-l,k+l)/k); 

T[k]  =  T[k-l]  +  t[k]; 
k  =  k  +  1; 

}  /*  end  while  */ 

a  =  levin(); 

}  /*  end  main  */ 

double  levinO 

{ 

double  P[MAX+1],  Q[MAX+1]; 
double  commonTerm; 
intN; 

P[0]  =  Q[0]  =  0.0; 

for  (N  =  1 ;  N  <=  MAX;  N++) 
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{ 

for  (k  =  1;  k  <=  N;  k++) 

{ 

cotnmonTerm  =  pow(-l,k)*pow(k,N-l)/tpc]*bico(N,k); 

P[k]  =  P[k-1]  +  conunonTenn*T[k]; 

Q[k]  =  Q[k-1]  +  commonTerm; 

} 

printf("\n%d  %32.3  Ilf ",  k-1,  log(2)  -  fabs(log(2)  -  P[N]/Q[N])); 

}  /*  end  for  N  loop  */ 

printf("\n\n"); 

return  (P[MAX]/Q(MAX]); 

} 

double  bico(int  n,  int  k) 

{  double  factln(int  n); 

return  floor(0.5+ejfp(factln(n)-factln(k)-factln(n-k))); 

} 

double  factln(int  n) 

{  double  ganunln(float  xx); 

static  float  a[101]; 
if  (n  <=  1)  return  0.0; 

if  (n  <=  100)  return  a[n]  ?  a[n] :  (a[n]=gainniln(n+1.0)); 
else  return  ganunln(n+1.0); 

} 

double  ganunln(float  xx) 

{  double  x,y,tmp,ser, 

static  double  cof[6]={76.18009172947146,-86.50532032941677, 
24.01409824083091,-1.231739572450155, 
0.1208650973866179e-2,-0.5395239384953e-5}; 

intj; 

y=x=xx; 

tnip=x+5.5; 

tmp  —  (x-K).5)*log(tmp); 

ser=  1.000000000 1900 15; 

for  0=0;j<=5  J-H-)  ser  +=  coflj]/++y; 

return  -tmp+-log(2.50662827463 10005*ser/x); 

} 
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APPENDIX  F.  MATLAB  CODE  TO  PRODUCE  PLOT  OF  CONVERGENCE  OF 

ALTERNATING  HARMONIC  SERIES 


/*  Thesis  Program  to  produce  data  for  plot  of  Levin  method  on  the  convergence  of  alternating  */ 
/*  harmonic  series  to  ln(2)  */ 

/*  LT  Brian  Johnson  */ 

/*  Software:  Matlab  Ver.  4.2C  */ 

/*  File  Name:  ln2comp.m  */ 

load  ln2graf  -ascii 
load  ln21evgr  -ascii 
X  =  ln2graf(:,l); 
y=ln2grafl;:,2); 

x2  =  ln21evgr(:,l); 
y2  =  ln21evgr(:,2); 

B  =  ones(l,71); 
figure(l) 

C  =  0.99*B; 

plot(x-l,y/log(2),x-l,C;:',  x2-l,y2/log(2)) 

xlabel('Number  of  iterations') 

ylabel('Partial  Sum  /  ln2') 

title  ('Levin  Method  vs  Infinite  Sum  for  ln(2)') 

axis([0  70  0  1]) 

gtext('0.99  — >') 
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APPENDIX  G.  ANSI  C  CODE  TO  PRODUCE  DATA  FOR  EXTINCTION 
EFFICIENCY  FACTOR  AND  CONVERGENCE  RATE  VS  SIZE  PARAMETER 


/*  Thesis  Program  to  produce  data  for  plotting  vs  size  parameter  */ 
/*  LT  Brian  Johnson  */ 

/*  Compiler:  Borland  C++  Ver.  5.0  */ 

/*  File  Name:  Miegraph.c  */ 

#include  <stdio.h> 

#include  <stdlib.h> 

#include  <stddef.h> 

/*  #mclude  <iostream.h>  */ 

#include  <math.h> 

#include  "nrutil.h" 

#include  "nrutil.c" 

#defineNRANSI 
#defineEPS  l.Oe-10 
#defineFPMIN  l.Oe-30 
#defineMAXIT  10000 
#defineXMIN2.0 
#define  PI  3.141592653589793 
#defme  RTPI02  1 .2533 141 
#defmeNUSEl  5 
#defmeNUSE2  5 

struct  bessVec  { 
double  rj,ry,rjp,  ryp; 

}; 


struct  complex  { 
double  real,  imag; 

}; 

void  assignBessels(struct  bessVec*  besselsPtr,  double  rjdummy,  double  rydummy,  double  rjpdummy, 
double  rypdummy) 

{ 

besselsPtr  ->  rj  =  rjdummy; 
besselsPtr  ->  ry  =  rydummy; 
besselsPtr  ->  rjp  =  rjpdummy; 
besselsPtr  ->  ryp  =  rypdummy; 

} 


struct  bessVec  sphbes(double,  int); 
struct  bessVec  bessjy(double,  double); 

struct  bessVec  sphBessels  n; 
struct  bessVec  bessels  nPlusHalf; 
struct  bessVec  bessels  nMinusHalf; 
struct  complex  a  n,  b  n; 
struct  complex  zeta  n; 
struct  complex  zeta  nMinusl; 
struct  complex  denom; 


/*  Function  Prototype  */ 

/*  Function  Prototype  */ 

/*  j's,  y s  &  derives  for  Ricatti-Bes  */ 
/*  for  calculating  J_n  +/-  half*/ 

/*  Mie  coefficients  */ 

/*  Ricatti-Bessel  fii  */ 

/*  denominator  in  Mie  coefficients  */ 
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double  x;  /* 

double  m;  /* 

double  J_nMmusHalf,  J  nPlusHalf;  /* 

double  A_factor;  /* 

double  psi_nMinusl;  /* 

double  psi_n; 

double  chi_n;  /* 

double  nOverx; 
double  aBraces,  bBraces; 
double  numer;  /* 

double  realDenom;  /* 

double  sumQextOld,  suni_Qext;  /* 

double  CLext[900][6];  /* 

mtnum[900][6];  /* 

double  conv  =  l.Oe-8; 


size  parameter  */ 

refractive  index  (non-absoibing  spheres)  */ 
for  calculating  A  factor  */ 
psi-prime  divided  1^  psi  */ 

Ricatti-Bessel  fri  1st  kind,  initialized  */ 


numerator  in  Mie  coefficients  */ 
rationalized  denominator  *! 

Extinction  efficiency  factor  */ 
array  of  Q_ext's,  fii's  of  m  and  x  */ 
number  of  iterations  to  converge,  f(x,m)  */ 


Ricatti-Bessel  fh  2nd  kind  */ 


inti; 

intj; 

intn; 


/*  counter  for  x  loop  */ 

/♦  counter  for  mArray  loop  */ 

/*  #  iterations  to  converge,  also  order  of  bessel  function  */ 


main() 

{ 

double  mArray[6]  =  {1.25,  1.33, 1.44,  1.50,  1.55,  2.0}; 


for(j  =  0;  j  <=  5;  j++)  { 


m  =  mArraylj]; 


for(x  ==  0.01;  x  <=9.0;  x  +=  0.01)  { 

n  =  0;  /*  initialize  n  for  each  x  */ 

bessels_nPlusHalf  =  bessjy(m*x,  0.5);  /♦  initialize  for  A  factor  */ 

JnPlusHalf  =  besselsnPlusHalf.ij; 

psi_n  =  sin(x); 

zeta_n.real  =  psi_n;  /*  initializing  zeta  */ 

zeta_n.imag  =  cos(x); 

sum_Qext=  0.0; 


do  { 

n  =  n  +  1; 

J_nMinusHalf  =  J_nPlusHalf; 
psi^nMinusl  =  psin; 
zeta_nMinusl  =  zeta_n; 
sumQextOld  =  sum_Qext; 

bessels_nPlusHalf  =  bessjy(m*x,  n  +  0.5);  /*  calculating  A  factor  */ 

J_nPlusHalf  =  bessels_nPlusHalf.ij; 

A_factor  =  (J_nMinusHalPJ_nPlusHalf)  -  (n/(m*x)); 
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/*  j's  and  y's  for  Ricatti  Bess  fii's  */ 


sphBesselsn  =  sphbes(x,  n); 
psi_n  =  X  *  sphBessels_n.q; 
chi  n  =  -X  *  sphBessels_n.iy; 
zeta_n.real  =  psin; 
zeta_n.imag  =  chi_n; 

nOverx  =  n/x; 

aBraces  =  (A_fector  /  m  +  nOverx); 
bBraces  =  (A_factor  *  m  +  nOverx); 

numer  =  aBraces*psi_n  -  psi_nMmusl; 

denom.real  =  (aBraces  *  zeta_n.real)  -  zetanMinusl.real; 
denom.iinag  =  (aBraces  *  zeta_n.unag)  -  zeta_nMinusl.iniag; 
realDenom  =  (denom.real’^denom.real  +  denom.imag’'‘denom.itnag); 

a_n.real  =  (numer  *  denom.real)/realDenom; 
a_n.imag  =  (numer  *  denom.imag)/realI>enom; 

numer  =  bBraces*psi_n  -  psinMnusl ; 

denom.real  =  (bBraces  *  zeta_n.real)  -  zeta_nMinusl.real; 
denom.imag  =  (bBraces  *  zeta_n.imag)  -  zetanMinusl.imag; 
realDenom  =  (denom.real'''denom.real  +  denom.imag*denom.imag); 

b_n.teal  =  (numer  *  denom.real)/realDenom; 
b_n.imag  =  (numer  *  denom.imag)/reaIDenom; 

sumQext  +=  (2*n  +  l)*(a_n.real  +  b_n.real); 

)  while  (fabs(sum_(Jext  -  sum_<3ext01d)  >  conv);  /*  end  do-while  loop  */ 

Q_ext[(int)(x*  100-1)]  [j]  =  2.0/(x*x)*sum_(Jext; 
mun[(mt)(x*  100-1)]  [j]  =  n; 

}  /*  end  iimer  x-loop  */ 

}  /*  end  outer  j-loop  (for  various  m  values)  */ 

for  (i  =  0;  i  <  899;  i++)  { 
printf("\n%f  ",  (i+1.0)/100.0); 


for(j  =  0;j<=5;j++) 
printf("%f  ",  CLext[i]D]); 

/*  printf("%d  ",  num[i]D]);  */ 


}  /*  end  i-j  loop  */ 


return  0; 

} 

struct  bessVec  sphbes(double  x,  int  n) 
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struct  bessVec  tempsphbessels; 
struct  bessVec  bessels; 

double  sj,  sy,  sjp,  syp; 
void  nrerror(char  error_text[]); 
double  factor,  order; 

if  (n  <  0 II X  <=  0.0)  urerrorC'bad  arguments  in  sphbes"); 
order=n+0.5; 

bessels  =  bessjy(x,order); 

factor=RTPI02/sqrt(x); 

sj=factor*bessels.ij; 

sy=factor*bessels.iy; 

sjp=factor*bessels.5p-(sj)/(2.0*x); 

syp=factor*bessels.ryp-(sy)/(2.0*x); 

assignBessels(&tempsphbessels,  sj,  sy,  sjp,  syp); 

return  tempsphbessels; 


} 

struct  bessVec  bessjy(double  x,  double  xnu) 

{ 

void  nrerror(char  error_text[]); 
struct  bessVec  tempBessels; 
double  rj,  ry,  gp,  lyp; 

voidbeschb(double  x,  double  *gaml,  double  *gam2,  double  *gampl, 
double  *gammi); 
int  i,isign,l,nl; 

double  a,b,br,bi,c,cr,ci,d,del,dell,den,di,dlr,dli,dr,e,f,fect,fact2, 

facts  ,flf,gam,gaml,gam2,gammi,gampl,h,p,pimu,pimu2,q,r,ijl, 

ijll,ijmu,ijpl,ijpl,qtemp,iyl,rymu,iymup,rytemp,sum,suml, 

temp,w,x2,xi,xi2,xmu,xmu2; 


if  (x  <=  0.0  II  xnu  <  0.0)  nrerror("bad  arguments  in  bessjy"); 

nl=(x  <  XMIN  ?  (int)(xnu+0.5)  :  IMAX(0,(int)(xnu-x+1.5))); 

xmu=xnu-nl; 

xmu2=xmu*xmu; 

xi=1.0/x; 

xi2=2.0*xi; 

w=xi2/PI; 

isign=l; 

h=xnu’'‘xi; 

if  (h  <  FPMIN)  h=FPMIN; 

b=xi2*xnu; 

d=0.0; 

c=h; 
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for(i=l;i<=MAXIT;i-H-)  { 
b  +=  xi2; 
d=b-d; 

if  (fabs(d)  <  FPMIN)  d=FPMIN; 
c=b-1.0/c; 

if  (fabs(c)  <  FPMIN)  c=FPM[N; 

d=1.0/d; 

del=c*d; 

h=del*h; 

if  (d  <  0.0)  isign  =  -isign; 
if  (fabs(del-1.0)  <  EPS)  break; 

} 

if  (i  >  MAXTT)  nrerror("x  too  large  in  bessjy;  tiy  asymptotic  expansion"); 

rjl=isign*FPMIN; 

qpl=h*ijl; 

rjpl=rjpl; 

fact=xnu*xi; 

for  (l=nl;l>=l;l-)  { 

qtemp=fact*ijl+qpl; 
fact  -=  xi; 

ijpl=fact*ijtemp-ijl; 

ijl=ijtemp; 

} 

if(ql==0.0)ql=EPS; 

f=ijpl/ijl; 

if(x<XMIN){ 

x2=0.5*x; 

pimu=PI*xmu; 

fact  =  (fabs(pimu)  <  EPS  ?  1.0  ;  pimt]/sin(punu)); 

d  =  -log(x2); 

e=xmu*d; 

fact2  =  (fabs(e)  <  EPS  ?  1.0  :  sinli(e)/e); 

beschb(xmu,&gaml,&gani2,&^mpl,&ganuni); 

fif=2.0/PI*fact*(gaml  *cosh(e)+gam2*fact2*d); 

e=exp(e); 

p=e/(gampl*PI); 

q=1.0/(e*PI*gamnii); 

pimu2==0.5*pimu; 

fact3  =  (fabs(pimu2)  <  EPS  ?  1.0  :  sin(pimu2)/pimu2); 

r=PI*pimu2*fact3  *fact3 ; 

c=1.0; 

d=-x2*x2; 

sum=fF+r*q; 

suml=p; 

for(i=l;i<=MAXIT;i++)  { 

fiP=(i*fiF+p+q)/(i*i-xmu2); 

c  *=  (d/i); 

p  /=  (i-xmu); 

q  /=  (i+xmu); 

del=c*(ff+r*q); 

sum  +=  del; 

dell=c*p-i*del; 
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suml  +=  dell; 

if  (fabs(del)  <  (1.0+fabs(sum))*EPS)  break; 


} 

if  (i  >  MAXIT)  nrerror('T5essy  series  failed  to  converge"); 

rymu  =  -sum; 

ryl  =  -siunl*xi2; 

rymup=xmu*xi*rymu-ry  1 ; 

ijmu=w/(iymup-f''rymu); 

}  else  { 

a=0.25-xmu2; 
p  =  -0.5*xi; 
q=1.0; 
br=2.0*x; 
bi=2.0; 

fact=a*xi/(p*p+q*q); 

cr=br+q*fact; 

ci=bi+p*fact; 

den=br*br+bi*bi; 

dr=br/den; 

di  =  -bi/den; 

dlr=cr*dr-ci*di; 

dli=cr*di+ci*dr; 

temp=p*dlr-q*dli; 

q=p*dli+q*dlr; 

p=temp; 

for  (i=2;i<=MAXIT;i-H-)  { 
a+=2*(i-l); 
bi  +=  2.0; 
dr=a*dr+br; 
di=a*di+bi; 

if  (fabs(dr)+M)s(di)  <  FPMIN)  dr=FPMIN; 

fact=a/(cr*cr+ci*ci); 

cr=br+cr*fact; 

ci=bi-ci*fact; 

if  (fabs(cr)+fabs(ci)  <  FPMEN)  cr=FPMIN; 

den=dr*dr+di*di; 

dr/=den; 

di  /=  -den; 

dlr=cr*dr-ci*di; 

dli=cr*di+ci*dr; 

temp=p*dlr-q*dli; 

q=p*dli+q*dlr; 

p=temp; 

if  (fabs(dlr-1.0)+fabs(dli)  <  EPS)  break; 

} 

if  (i  >  MAXIT)  nrerror("cf2  failed  in  bessjy"); 
gam=(p-f)/q; 

ijmu=sqrt(w/((p-f)*gam+q)); 

ijmu=SIGN(ijmu,ijl); 

rymu=ijmu*gam; 

rymup=tymu*(p+q/pm); 

iyl=xmu*xi*iymu-rymup; 


} 
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fact=ijmu/ql; 
q=rjll*fact; 
ijp=ijpl*fact; 
for  (i=l;i<=iil;i-H-)  { 

iytemp=(xmu+i)*xi2*iyl-rymu; 

rymu=iyl; 

ryl=iytemp; 

} 

iy=iymu; 

ryp=xnu*xi*iymu-ry  1 ; 
assignBessels(&tempBessels,  ij,  ry,  ijp,  lyp); 
return  tempBessels; 

} 

voidbeschb(double  x,  doable  *gaml,  double  *gain2,  double  *gampl,  double  *gamnu) 

{ 

double  chebev(double  a,  double  b,  double  c[],  int  m,  double  x); 
double  xx; 

static  double  cl[]  =  { 

-1. 142022680371 172e0,6.51651 1267076e-3, 
3.08709017308e-4,-3.470626964e-6,6.943764e-9, 
3.6780e-ll,-1.36e-13}; 
static  double  c2[]  =  { 

1.843740587300906e0,-0.076852840844786e0, 

1.271927136655e-3,-4.971736704e-6,-3.3126120e-8, 

2.42310e-10,-1.70e-13,-1.0e-15>; 

xx=8.0*x*x-1.0; 

*gaml=chd)ev(-1.0,1.0,cl,NUSEl,xx); 

*gain2=clid)ev(-1.0,1.0,c2,NUSE2,xx); 

*gampl=  *gain2-x*(*gaml); 

*gammi=  ♦gam2+x*(*gaml); 

} 

double  chebev(double  a,  double  b,  double  cQ,  int  tn,  double  x) 

{ 

void  nrerroifchar  errortextQ); 
double  d  =  0.0,  dd  =  0.0,  sv,  y,  y2; 
intj; 

if  ((x-a)*(x-b)  >  0.0)  nrerror("x  not  in  range  in  routine  chdiev''); 
y2  =  2.0*(y=(2.0*x-a-b)/(b-a)); 
for(j  =  m-l;j>=  l;j-){ 
sv  =  ^ 

d  =  y2*d-dd  +  c|j]; 
dd  =  sv; 

} 

return  y*d  -  dd  +  0.5*c[0]; 

} 


#undefEPS 
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#undefFPMIN 

#undefMAXIT 

#undefXMIN 

#undefPI 

#undefNRANSI 

#imdefNUSEl 

#undefNUSE2 

#imdefRTPI02 


APPENDIX  H.  MATLAB  CODE  TO  PLOT  OUTPUT  OF  CODE  LISTED 

IN  APPENDIX  G. 


load  Qout  -ascii 
X  =  Qout(:,l); 
yl  =  Qout(:,2); 
y2  =  Qout(:,3); 
y3  =  Qout(:,4); 
y4  =  Qout(:,5); 
y5  =  Qout(:,6); 
y6  =  Qout(:,7); 

step  =  0.25; 
xinc  =  step:step;8.75; 
for  i  =  1  :(length(xmc)-l) 
ylinc(i)  =  yl(100*xinc(i)+l); 
y2inc(i)  =  y2(100*xinc(i)+l); 
y3inc(i)  =  y3(100*xinc(i)+l); 
y4inc(i)  =  y4(100*xinc(i)+l); 
y5inc(i)  =  y5(100*xinc(i)+l); 
y6inc(i)  =  y6(100*xinc(i)+l); 
end 

plot(x,yl,x,y2,x,y3,x,y4,x,y5,x,y6,  xinc,ylinc,'.',  xinc,y2inc,’x',  xinc,  y3inc,'*',  xinc,y4inc,'+', 

xinc,y5inc,'o',  xinc,  y6inc,'.') 

axis([0,9,0,6]) 

xlabel('Size  parameter  x'),  ylabel('Qext') 

gtext('m  =  2') 

gtext('m  =  1.55') 

gtext('m  =  1.50') 

gtext('m  =  1.44') 

gtext('m  =  1.33') 

gtext('m  =  1.25') 

load  Qoutn  -ascii 
X  =  Qoutn(:,l); 
yl  =  Qoutn(:,2); 
y2  =  Qoutn(:,3); 
y3  =  Qoutn(:,4); 
y4  =  Qoutn(:,5); 
y5  =  Qoutn(:,6); 
y6  =  Qoutn(;,7); 

step  =  0.25; 

xinc  =  step:step:8.75; 

for  i  =  l:(length(xinc)-l) 
ylinc(i)  =  yl(100*xinc(i)+l); 
y2inc(i)  =  y2(100’*‘xinc(i)+l); 
y3inc(i)  =  y3(100'''xinc(i)+l); 
y4inc(i)  =  y4(100’*'xinc(i)+l); 
y5inc(i)  =  y5(100'''xinc(i)+l); 
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y6mc(i)  =  y6(100*xinc(i)+l); 
end 

plot(x,yl,'k',x,y2,'m',x,y3,’c',x,y4,’r',x,y5,'g’,x,y6,'b’,  xinc,ylinc,'.k',  xinc,y2inc,’xm',  xinc,  y3inc,'*c’, 
xinc,y4inc,'+r',  xinc,y5inc,'og',  xinc,  y6inc,'.b') 

xIabel('Size  parameter  x'),  ylabel('#  iterations  to  converge’) 
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APPENDIX  1.  ANSI  C  CODE  FOR  CALCULATING  EXTINCTION 

EFFICIENCY  FACTOR 


/*  Thesis  Program  for  calculating  as  a  fimction  of  user-defined  size  parameter  */ 
/*  LT  Brian  Johnson  */ 

/*  Compiler:  Borland  C++  Ver.  5.0  */ 

/*  File  Name:  mie.c  */ 


#include  <stdio.h> 

#include  <stdlib.h> 

#include  <stddef.h> 

#include  <math.h> 

#include  "nrutiLh" 

#include  "nrutilc" 

#defineNRANSI 

#define  EPS  l.Oe-10 

#defineFPMIN  l.Oe-30 

#define  MAXIT  10000 

#defineXMIN2.0 

#define  PI  3.141592653589793 

#defme  RTPI02  1.2533141 

#defineNUSEl  5 

#defineNUSE2  5 

struct  bessVec  { 
double  rj,  ry,  rjp,  ryp; 

}; 


Struct  complex  { 
double  real,  imag; 

}; 


void  assignBessels(struct  bessVec*  besselsPtr,  double  rjdummy,  double  rydummy,  double  rjpdummy, 
double  rypdummy) 

{ 

besselsPtr  ->  rj  =  rjdummy; 
besselsPtr  ->  ry  =  rydummy; 
besselsPtr  ->  rjp  =  rjpdummy; 
besselsPtr  ->  ryp  =  rypdummy, 

} 

struct  bessVec  sphbes(double,  int);  /*  Function  Prototype  ♦/ 

struct  bessVec  bessjy(double,  double);  /*  Function  Prototype  */ 

struct  bessVec  sphBessels  n; 

struct  bessVec  bessels  nPlusHalf; 

struct  bessVec  bessels  nMinusHalf; 

struct  complex  a  n,  b  n;  /*  Mie  coefBcients  */ 

struct  complex  zeta  n; 

struct  complex  zeta  nMinusl ; 

struct  complex  denom;  /*  denominator  in  Mie  coefficients  */ 
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double  X  =  109.0;  /*  size  parameter  */ 

double  m  =  2.0;  /*  refiactive  index  (non-absorbing  spheres)  */ 

double  J_nMinusHalf,  JnPlusHalf; 

double  A_factor; 

double  psi_nMinusl;  /*  Ricatti-Bessel  fii  1st  kind,  initialized  */ 

double  psi_n; 
double  chi_n; 
double  nOverx; 
double  aBraces,  bBraces; 
double  numer; 
double  realDenom; 

double  sum  QextOld,  sum  Qext,  Q_ext; 
double  conv  =  l.Oe-8; 


/♦  Ricatti-Bessel  fh  2nd  kind  */ 


/*  numerator  in  Mie  coefficients  */ 

/*  rationalized  denominator  */ 

/*  Extinction  efficiency  factor  */ 


int  n  =  0;  /*  loop  counter,  also  order  of  bessel  hmction  */ 


main() 

{ 

bessels^nPlusHalf = bessjy(m*x,  0.5);  /*  initialize  for  A  factor  */ 

J_nPlusHalf  =  bessels_nPlusHalf.ij; 

psi_n  =  sin(x); 

zeta_n.real  =  psi_n;  /*  initializing  zeta  */ 

zeta_n.imag  =  cos(x); 

sum_Qext=  0.0; 

do  { 

n  =  n+  1; 

J_iiMmusHaIf  =  J_nPlusHalf; 
psinMinusl  =  psin; 
zeta_nMinusl  =  zeta_n; 
sumQextOld  =  sumQext; 

bessels  nPlusHalf  =  bessjy(m*x,  n  +  0.5);  /*  calculating  A  factor  */ 

J_iiPlusHalf  =  bessels_nPlusHalf.q; 

A_factor  =  (JnMinusHaffiJnPlusHalf)  -  (n/(m*x)); 

sphBessels_n  ==  sphbes(x,  n);  /*  j's  and  y's  for  Ricatti  Bess  fit's  */ 

psi_n  =  X  *  sphBessels_n.rj; 

chin  =  -X  ♦  sphBessels_n.ry; 

zeta_n.real  =  psi_n; 

zeta_n.tmag  =  chin; 

nOverx  =  n/x; 

aBraces  =  (A_factor  /  m  +  nOverx); 
bBraces  =  (A_factor  *  m  +  nOverx); 

numer  =  aBraces*psi_n  -  psi_nMinusl; 

denom.real  =  (aBraces  *  zeta_n.real)  -  zeta_nMinusl.real; 
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denom.iinag  =  (aBraces  *  zeta_n.iinag)  -  zetanMinusl.imag; 
realDenom  =  (denom.real*denom.reaI  +  denom.iinag*^denom.iiiiag); 

an.real  =  (numer  *  denom.real)/realDenom; 
a_n.imag  =  (numer  *  denom.imag)/reaIDenom; 

numer  =  bBraces*psi_n  -  psinMinusl; 

denom.real  =  (bBraces  *  zeta_n.real)  -  zeta_nMinusl.real; 
denom.imag  =  (bBraces  *  zeta_n.imag)  -  zeta  nMinusl.imag; 
realDenom  =  (denom.real*denom.real  +  denom.imag*denom.imag); 

b_n.real  =  (numer  *  denom.real)/realDenom; 
b_n.imag  =  (numer  *  denom.imag)/realDenom; 

sumQext  +=  (2*n  +  l)*(a_n.real  +  b_n.real); 

}  while  (fabs(sum_Qext  -  sum_Qext01d)  >  conv); 

Q^ext  =  2.0/(x*x)*sum_CJext; 

printf("\n\nThe  Extinction  EflBciency  Factor  is:  Qext  =  %r,  Q_ext,  "\n\n"); 
printf("\nNumber  of  iterations:  %d\n",  n); 
printf(''\nx  =  %3.2f  m  =  %1.2f  \n\n",  x,  m); 
return  0; 

} 

struct  bessVec  sphbes(double  x,  int  n) 

{ 

struct  bessVec  tempsphbessels; 
struct  bessVec  bessels; 

double  sj,  sy,  sjp,  syp; 
void  nrerTor(char  error_textD); 
double  factor,  order; 

if  (n  <  0  II X  <=  0.0)  nrerror("bad  arguments  in  sphbes"); 
order=n+0.5; 

bessels  =  bessjy(x,order); 

fector=RTPI02/sqrt(x); 

sj=factor*bessels.ij; 

sy=factor*bessels.ry; 

sjp=factor*bessels.rjp-(sj)/(2.0*x); 

syp=factor*bessels.typ-(sy)/(2.0*x); 

assignBessels(&tempsphbessels,  sj,  sy,  sjp,  syp); 

return  tempsphbessels; 

} 

struct  bessVec  bessjy(double  x,  double  xnu) 

{ 
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void  iirerror(char  error_text[]); 
struct  bessVec  tempBessels; 
double  q,  ly,  qp,  ryp; 

void  beschb(double  x,  double  *gaml,  double  *gain2,  double  *gampl, 
double  *gamTni); 
int  i,isign,l,nl; 

double  a,b,br,bi,c,cr,ci,d,del,del  1  ,den,di,dlr,dli,dr,e,f,fact,fact2, 

facts  ,flF, gam, gaml,gam2,gainnii,gampl,h,p,pimu,pimu2,q,r,ql, 

qll,qmu,qpl,qpl,qtemp,iyl,rymu,iymup,rytemp,sum,suml, 

temp,w,x2,xi,xi2,xmu,xmu2; 


if  (x  <=  0.0  II  xnu  <  0.0)  nrenor("bad  arguments  in  bessjy"); 

nl=(x<XMIN  ?  (int)(xnu-K).5) ;  IMAX(0,(int)(xnu-x+1.5))); 

xmu=xnu-nl; 

xmu2=xmu*xmu; 

xi=1.0/x; 

xi2=2.0*xi; 

w=xi2/PI; 

isign=l; 

h=xnu*xi; 

if  (h  <  FPMIN)  h=FPMIN; 

b^xi2*xnu; 

d=0.0; 

c=h; 

for  (i=l;i<=MAXrr;i-H-)  { 
b+=  xi2; 

d=b-d; 

if  (fabs(d)  <  FPMIN)  d=FPMIN; 
c=b-1.0/c; 

if  (fabs(c)  <  FPMIN)  c=FPMIN; 

d=1.0/d; 

del=c*d; 

h=del*h; 

if  (d  <  0.0)  isign  =  -isign; 
if  (fabs(del-1.0)  <  EPS)  break; 

} 

if  (i  >  MAXTT)  nrerror("x  too  large  in  bessjy;  try  a^mptotic  expansion"); 

rjl=isign*FPMIN; 

qpl=h*ql; 

qii=rji; 

qpl=qpl; 

fact=xnu*xi; 

for  (l=nl;l>=l;l--)  { 

qtemp=fact*ql+qpl; 
fart  -=  xi; 

qpl=fart*qtemp-ql; 

ql=qtemp; 

} 

if(ql==0.0)  ql=EPS; 
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f=qpl/ijl; 

if(x<XMIN){ 

x2=0.5*x; 

pimu=PI*xmu; 

fact  =  (fab5(piinu)  <  EPS  ?  1,0  :  punu/sm(pimu)); 

d=-log(x2); 

e=xniu*d; 

fact2  =  (fabs(e)  <  EPS  ?  1.0  :  siiih(e)/e); 

beschb(xmu,&gaml,&gam2,&gainpl,&gammi); 

fif=2.0/Pl*fiict*(gaml*cosh(e)+gain2*fact2*d); 

e=exp(e); 

p=^(gampl*PI); 

q=  1 .0/(e*Pl*gammi); 

piinii2=0.5*pimu; 

facts  =  (f9bs(pimu2)  <  EPS  ?  1.0  :  sin(pimu2)/pimu2); 

r=PI*pimu2*fact3  *fect3 ; 

c=1.0; 

d  =  -x2*x2; 

sum=fF+r*q; 

suml=p; 

for  (i=l;i<=MAXIT;i-H-)  { 

ff=(i*ff+pfq)/(i*i-xmu2); 
c  *=  (d/i); 
p  /=  (i-xmu); 
q/=  (i+xmu); 

(fcl=c*(ff+r*q); 
sum  +=  del; 
dell=c*p-i*del; 
suml  +=  dell; 

if  (fabs(del)  <  (1.0+fabs(sum))*EPS)  bie^ 

} 

if  (i  >  MAXIT)  nrerrorCbessy  series  failed  to  converge"); 

rymu  =  -sum; 

lyl  =  -suml*xi2; 

rymiip=xmu*xi*rymu-ryl; 

ijmu=w/(rymup-f'‘iymu); 

}  else  { 

a=0.25-xmu2; 
p=  -0.5*xi; 
q=1.0; 
br=2.0*x; 
bi=2.0; 

fact=a*xi/(p*p+-q*q); 

cr=br+q*fact; 

ci=bi+p*fact; 

den=br*brfbi*bi; 

dr=br/den; 

di  =  -bi/den; 

dlf=cr*dr-ci*di; 

dli=cr*di+ci*dr; 

temp=p*dlr-q*dli; 

<f=p*dli+q*dlr; 

p=lemp; 
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for  (i=2;i<=MAXIT;i-H-)  { 
a+=2*(i-l); 
bi  +=  2.0; 
dr=a*di:+br; 
di=a*di+bi; 

if  (fabs(dr)+fabs(di)  <  FPMIN)  dr=FPMIN; 

fact=a/(cr*cr+ci*ci); 

cr=br+ca-*fact; 

ci=bi-ci*fact; 

if  (fabs(cr)+fabs(ci)  <  FPMIN)  cr=FPMrN; 

den=dr*dr+di*di; 

dr  /=  den; 

di  /=  -den; 

dlr=cr*dr-ci*di; 

dli=cr*di+ci*dr; 

tenip=p*dlr-q*dli; 

q=p*dli+q*dlr; 

p=temp; 

if  (fabs(dlr-1.0)+fabs(dli)  <  EPS)  break; 

> 

if  (i  >  MAXIT)  nrerror("cf2  failed  in  bessjy"); 
gam=(p-f)/q; 

ijmu=sqrt(w/((p-0*gam+q)); 

qmu=SIGN(ijmu,ql); 

rymu=rjmu*gam; 

iyinup=rymu*(p+q/gani); 

ryl=xmu*xi*rymu-iymup; 

} 

fact=qmn/ql; 
ij=ijll*fect; 
rjp=qpl*fact; 
for  (i=l;i<=nl;i-H-)  { 

rytemp^(xmu+i)*xi2*ryl-rymu; 

iyinu=iyl; 

iyl=rytemp; 

} 

ty=iymu; 

ryp=xnu*xi*rymu-ryl ; 


assignBessels(&tempBessels,  q,  ry,  qp,  typ); 
return  tempBessels; 

} 

void  beschb(double  x,  double  *gaml,  double  *gain2,  double  *gampl,  double  *ganuiii) 

{ 

double  chebev(double  a,  double  b,  double  cQ,  int  m,  double  x); 
double  xx; 

static  double  clfl  =  { 

-1 . 142022680371 172e0,6.5 165 1 1267076e-3, 
3.08709017308e-4,-3.470626964e-6,6.943764e-9, 
3.6780e-ll,-1.36e-13}; 
static  double  c2Q  =  { 
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1.843740587300906e0,-0.076852840844786e0, 

1.271927136655e-3,-4.971736704e-6,-3.3126120e-8. 

2.42310e-10,-1.70e-13,-1.0e-15}; 


xx=8.0*x*x-1.0; 

*gaml=chebev(-1.0,1.0,cl,NUSEl,xx); 

*gam2=chd)ev(-l  .0, 1 .0,c2,NUSE2,xx); 

*gampl=  *gam2-x*(*gaml); 

*gammi=  *gam2+x*(*gaml); 

> 

double  chebev(double  a,  double  b,  double  cQ,  int  m,  double  x) 

{ 

void  nrerror(char  errorJextQ); 
double  d  =  0.0,  dd  =  0.0,  sv,  y,  y2; 
intj; 

if  ((x-a)*(x-b)  >  0.0)  nreiTor("x  not  in  range  in  routine  chebev") 
y2  =  2.0*(y=(2.0*x-a-b)/(b-a)); 
forO  =  m-l;j>= 
sv  =  ^ 

d  =  y2*d-dd  +  c[j], 

dd  =  sv; 

} 

return  y*d-dd  +  0.5*c[0]; 

} 

#undefEPS 

#undefFPMIN 

#undefMAXIT 

#undefXMIN 

#undefPI 

#undefNRANSI 

#undefNUSEl 

#undefNUSE2 

#undefRTPI02 
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APPENDIX  J.  ANSI  C  CODE  FOR  LEVIN  ALGORITHM  OPERATING  ON 
EXTINCTION  EFFICIENCY  FACTOR 


/*  Thesis  Program  for  applying  the  Levin  method  on  the  extinction  efficiency  factor  and  determining  */ 
/*  the  resulting  rate  of  convergence  */ 

/*  LT  Brian  Johnson  */ 

/*  Compiler:  Borland  C++  Ver.  5.0  */ 

/*  File  Name:  mielevin.c  */ 

#include  <stdio.h> 

#include  <stdlib.h> 

#include  <stddef.h> 

#include  <math.h> 

#include  "nrutil.h” 

#include  ”nrutil.c” 

#defineNRANSI 
#defineEPS  l.Oe-10 
#defineFPMIN1.0e-30 
#defme  MAXIT  10000 
#define  XMIN  2.0 
#define  PI  3.141592653589793 
#define  RTPI02  1.2533141 
#defineNUSEl  5 
#defineNUSE2  5 
#define  MAX  1036 

struct  bessVec  { 
double  rj,  ry,  rjp,  ryp; 

}; 


struct  complex  { 
double  real,  imag; 

}; 


void  assignBessels(struct  bessVec*  besselsPtr,  double  rjdummy,  double  rydummy,  double  rjpdummy, 
double  rypdummy) 

{ 

besselsPtr  ->  rj  =  rjdummy; 
besselsPtr  ->  ry  =  rydummy; 
besselsPtr  ->  rjp  =  rjpdummy; 
besselsPtr  ->  ryp  =  rypdummy; 

} 


struct  bessVec  sphbes(double,  int); 
struct  bessVec  bessjy(double,  double); 

double  levinO; 
double  bico(int,  int); 
double  factln(mt); 
double  gammln(double); 


/*  Function  Prototype  */ 

/*  Function  Prototype  */ 

/*  Function  Prototype  */ 
/*  Fimction  Prototype 
I*  Function  Prototype  */ 
/*  Function  Prototype  */ 


struct  bessVec  sphBessels  n; 
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struct  bessVec  besselsnPlusHalf; 

struct  bessVec  bessels_iiMinusHalf; 

struct  complex  a_ii,  b_n;  /*  Mie  coefficients  */ 

struct  complex  zeta  n; 

struct  complex  zeta_nMinusl ; 

struct  complex  denom;  /*  denominator  in  Mie  coefficients  */ 


double  X  =  1000.0;  /*  size  parameter  *! 

double  m  =  2.0;  /*  refractive  index  (non-absorbing  spheres)  */ 

double  J  nMinusHalf,  J  nPlusHalf; 
double  A_factor; 
double  psi_nMinusl; 
double  psin; 
double  chi_n; 
double  nOverx; 
double  aBraces,  bBraces; 
double  numer; 
double  realDenom; 

double  sum  QextOld,  sum  Qext,  Q_ext; 
double  T[MAX+1],  t[MAX+l]; 
double  a; 


I*  Ricatti-Bessel  fh  1st  kind,  initialized  *! 


/*  Ricatti-Bessel  fii  2nd  kind  */ 


/*  munerator  in  Mie  coefficients  */ 
/*  rationalized  denominator  *! 

/♦  Extinction  efficiency  fector  */ 


int  n  =  0;  /*  loop  counter,  also  order  of  bessel  function  */ 

int  mArrayn  =  {1.25,  1.33, 1.44, 1.50,  1.55, 2.0}; 

intk; 

mainO 

{ 

T[0]=t[0]  =  0.0; 
n  =  0; 

printf("\n\ni  T_i  t_i  Pti]/Q[i]\n"); 

bessels_nPlusIfelf  =  bessjy(m*x,  0.5);  /*  initialize  for  A_factor  */ 

J_nPlusHalf  =  bessels_nPlusHalf.q; 

psi_n  =  sin(x); 

zeta_n.real  =  psi_n;  /*  initializing  zeta  ♦/ 

zeta_n.imag  =  cos(x); 

sum_Qext=  0.0; 

while  (n<=MAX) 

n  =  n+  1; 

J_nMinusHalf  =  J_nPlusHalf; 
psi_nMinusl  =  psi_n; 
zeta_nMinusl  =  zetan; 
sumQextOld  =  sumQext; 

bessels_nPlusHalf  =  bessjyf™**.  n  +  0-5);  I*  calculating  A  factor  */ 

JnPludlalf = besselsnPlusHalfq; 

A_&ctor  =  (J_nMinusHaffiJ_nPlusHalf)  -  (n/(m*x)); 
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sphBessels_n  =  sphbes(x,  n);  /*j's  and  y's  for  Ricatti  Bess  fh's*/ 

psi_n  =  X  *  sphBessels_n.q; 

chin  =  -X  *  sphBessels_n.ry; 

zeta_n.real  =  psi_n; 

zeta_n.imag  =  chi_n; 

nOverx  =  n/x; 

aBraces  =  (A_factor  /  m  +  nOverx); 
bBraces  =  (Afactor  *  m  +  nOverx); 

numer  =  aBiaces*psi_n  -  psi_nMinusl; 

denom.real  =  (aBraces  *  zeta_n.real)  -  zetanMinusl.real; 
denom.imag  =  (aBraces  *  zeta_n.iniag)  -  zeta_nMmusl.imag; 
realDenom  =  (denoni.real*denom.real  +  denom.iniag*denom.imag); 

a  n.real  =  (numer  *  denom.real)/realDenom; 
a_n.iniag  =  (numer  *  denom.imag)/realDenom; 

numer  =  bBraces*psi_n  -  psi_nMinusl; 

denom.real  =  (bBraces  *  zeta_n.real)  -  zeta_nMinusl.real; 
denom.imag  =  (bBraces  *  zeta_n.imag)  -  zetanMinusl.imag; 
realDenom  =  (denom.real*denom.real  +  denom.imag*denom.imag); 

b_n.real  =  (numer  *  denom.real)/realDenom; 
b_n.imag  =  (numer  *  denom.imag)/realDenom; 

t[n]  =  2.0/(x*x)*(2*n  +  l)*(a_n.real  +  b_n.real); 

T[n]  =  T[n-l]  +  t[n]; 

}  /*  end  while  loop  */ 

a  =  levinO; 

printf("\nx  =  %3.2f  m  =  %1.2f  \n\n",  x,  m); 


return  0; 

} 

struct  bessVec  sphbes(double  x,  int  n) 

{ 

struct  bessVec  tempsphbessels; 
struct  bessVec  bessels; 

double  sj,  sy,  sjp,  syp; 
void  nrerror(char  error_text[]); 
double  factor,  order; 

if  (n  <  0  II  x  <=  0.0)  nrerrorCTjad  arguments  in  sphbes"); 
order=n-K).5; 

bessels  =  bessjy(x,order); 
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factor=RTPI02/sqrt(x); 

sj=factor*bessels.ij; 

sy=factor*bessels.iy; 

sjp=factor*bessels.rji>(sj)/(2.0*x); 

syp=factor*bessels.ryp-(sy)/(2.0*x); 

assignBessels(&tempsphbessels,  sj,  sy,  sjp,  syp); 

return  tempsphbessels; 

} 

struct  bessVec  bessjy(double  x,  double  xnu) 

{ 

void  nrerror(char  enor_text[]); 
struct  bessVec  tempBessels; 
double  Tj,  ly,  ijp,  lyp; 

void  beschb(double  x,  double  *gaml,  double  *gam2,  double  *gampl, 
double  *ganimi); 
int  i,isign,l,nl; 

double  a,b,br,bi,c,cr,ci,d,del,dell,den,di,dlr,dli,dr,e,f, fact, fact2, 

fact3,fF,gam,gaml,gam2,ganuni,gampl,h,p,pimu,pimu2,q,r,ijl, 

ijll,rjmu,rjpl,qpl,qtemp,iyl,rymu,Tymup,rytemp,sum,siunl, 

temp,w,x2,xi,xi2,xmu,xmu2; 


if  (x  <=  0.0  II  xnu  <  0.0)  nrerror("bad  arguments  in  bessjy"); 

nl=(x  <  XMIN  ?  (int)(xnu+0.5)  ;  lMAX(0,(int)(xnu-x+1.5))); 

xmu=xnu-nl; 

xmu2=xmu*xmu; 

xi=1.0/x; 

xi2=2.0*xi; 

w=xi2/PI; 

isign=l; 

h=xnu*xi; 

if  (h  <  FPMIN)  h=FPMIN; 
b=xi2*xnu; 

<^0.0; 

c=h; 

for  (i=l;i<=MAXrr;i-H-)  { 
b+=  xi2; 
d=b-d; 

if  (fabs(d)  <  FPMIN)  d=FPMIN; 
c=b-1.0/c; 

if  (fabs(c)  <  FPMIN)  c=FPMIN; 

d=1.0/d; 

del=c*d; 

h=del*h; 

if  (d  <  0.0)  isign  =  -isign; 
if  (fabs(del-1.0)  <  EPS)  break; 


} 
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if  (i  >  MAXrr)  iireiTor("x  too  large  in  bessjy;  try  asymptotic  expansion"); 

ql=isign*FPMIN; 

ijpl=h*rjl; 

r(ll=rjl; 

ijpl=qpl; 

fact=xnu*xi; 

for  (l=nl;l>=l;l~)  { 

ijtemp=fact*ijl+ijpl; 
fact  -=  xi; 

gpl=fact*ijtemp-ql; 

ql=ijtemp; 

} 

if  (rjl  =  0.0)  ijl=EPS; 

f=rjpl/gl; 

if(x<XMIN){ 

x2=0.5*x; 

pimu=PI*xmu; 

fact  =  (fabs(pimu)  <  EPS  ?  1.0  :  pimii/sin(punu)); 

d=-log(x2); 

e=xmu*d; 

fect2  =  (fabs(e)  <  EPS  ?  1.0  :  sinh(e)/e); 

beschb(xmu,&gaml,&gam2,&gampl,&ganuni); 

fi=2.0/PI*fact*(gaml*cosh(e)+gain2*fact2*d); 

e=exp(e); 

pF=e/(gampl*PI); 

q=  1 .0/(e*PI*gamnu); 

pimu2=0.5*pimn; 

facts  =  (f^(punu2)  <  EPS  ?  1.0  :  sin(pimu2)/pimu2); 

r=Pl*pimu2*fact3  *fect3 ; 

c=1.0; 

d=-x2*x2; 

sum=fF+r*q; 

suml=p; 

for  (i=l;i<=MAXIT;i++)  { 

fP=(i*fiF+iH-q)/(i*i-xmu2); 
c*=(d/i); 
p  /=  (i-xmu); 
q  /=  (i+xmu); 
del=c*(ff+r*q); 
sum  +=  del; 
dell=c*p-i*del; 
suml  +=  dell; 

if  (fabs(del)  <  (1.0+febs(sum))*EPS)  break; 

} 

if  (i  >  MAXIT)  nrerror("bessy  series  failed  to  converge"); 

rymu  =  -sum; 

ryl  =  -suml*xi2; 

rymup=xmu*xi*rymu-ry  1 ; 

rjmu=w/(rymup-f'‘rymu); 

}  else  { 

a=0.25-xmu2; 

p=-0.5*xi; 

q=1.0; 
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br=2.0*x; 

bi=2.0; 

fact=a*xi/(p*p+q*q); 

cr=bH-q*fact; 

ci=bi+p*fact; 

den=br*br+bi*bi; 

dr=br/den; 

di  =  -bi/den; 

dlr=cr*dr-ci*di; 

dli=cr*di+d*dr, 

temp=p*dlr-q*dli; 

q=p*dli+q*dlr; 

p=temp; 

for  (i=2;i<=MAXIT;i++)  { 
a+=2*(i-l); 
bi  +=  2.0; 
dr=a*dr+-br; 
di=a*di+bi; 

if  (fabs(dr)+fabs(di)  <  FPMIN)  dr=FPMIN; 

fact=a/(cr*cr+d*ci); 

cr=br+cr*fiict; 

ci=bi-ci*fact; 

if  (fabs(cr)+febs(ci)  <  FPMIN)  cr=FPMIN; 

den=dr*dr+di*di; 

dr  /=  den; 

di  /=  -den; 

dlr=cr*dr-ci*di; 

dli=cr*di+ci*dr; 

temp=p*dlr-q*dli; 

q=p*dli+q*dlr, 

p^mp; 

if  (fabs(dlr-1.0)+fabs(dli)  <  EPS)  break; 

} 

if  (i  >  MAXTT)  nrerror("cf2  failed  in  bessjy"); 
gam=(p-f)/q; 

gmu=sqrt(w/((p-f)*gani+q)); 

gmu=SIGN(ijnni,ql); 

rymu=ijmu*gam; 

rymiqjF=iymu*(p+q/gam); 

ryl=xmu*xi*iymu-ryniiip; 

} 

fact=ijmii/ijl; 
q=qll*fact; 
qp=ijpl*fact; 
for  (i=l;i<=nl;i-H-)  { 

rytemp=(xmu+i)*xi2*ry  1  -rymu; 

iyinu=ryl; 

iyl=rytemp; 

} 

ry=iyniu; 

ryp=xnu*xi*iymu-iyl; 

assignBessels(&tempBessels,  ij,  ry,  rjp,  ryp); 
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return  tempBessels; 

} 

void  beschb(double  x,  double  *gaml,  double  *gain2,  double  *gampl,  double  *gainim) 

{ 

double  chebev(double  a,  double  b,  double  cQ,  int  m,  double  x); 
double  xx; 

static  double  cl[]  =  { 

-1. 142022680371 172e0, 6.5 1651 1267076e-3, 
3.08709017308e-4,-3.470626964e-6,6.943764e-9, 
3.6780e-ll,-1.36e-13}; 
static  double  c2[]  =  { 

1.843740587300906e0,-0.076852840844786e0, 

1.271927136655e-3,-4.971736704e-6,-3.3126120e-8, 

2.42310e-10,-1.70e-13,-1.0e-15}; 

xx=8.0*x*x-1.0; 

*gaml=chebev(-1.0,1.0,cl,NUSEl,xx); 

*gain2=chd)ev(-1.0,1.0,c2,NUSE2,xx); 

*gainpl=  *gain2-x*(*gaml); 

*ganuni=  *gani2+x*(*gaml); 


double  chebev(double  a,  double  b,  double  cQ,  int  m,  double  x) 

{ 

void  nrerror(char  error_textD); 
double  d  =  0.0,  dd  =  0.0,  sv,  y,  y2; 
intj; 

if  ((x-a)*(x-b)  >  0.0)  merror("x  not  in  range  in  routine  chebev"); 
y2  =  2.0*(y=(2.0*x-a-b)/(b-a)); 
forO  =  m-l;j>=l;j-){ 
sv  =  d; 

d  =  y2*d-dd  +  c01; 
dd  =  sv; 

} 

return  y*d  -  dd  +  0.5*c[0]; 


double  levinO 

{ 

double  P[MAX+1],  Q[MAX+1]; 
double  commonTerm; 
intN; 

PI0]  =  Q[0]  =  0.0; 

for  (N  =  1;  N  <=  MAX;  N++) 

{ 

for  Gc  =  1;  k  <=  N;  k-H-) 

{ 

conunonTerm  =  pow(-l,k)*pow(k,N-l)/t|k]*bico(N,k); 
PPc]  =  P[k-1]  +  coinmonTerm*T[k]; 
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Q[k]  =  Q[k-1]  +  commonTenn; 

} 

printf("\n%d  %f  %f  %f ",  k- 1,  T[N],  t[N],  P[N]/Q[N1); 

} 

printf("\n\n"); 

return  (P[MAX]/Q[MAX]); 

} 

double  bico(int  n,  int  k) 

{  double  factln(int  n); 

return  floor(0.5+exp{factln(n)-factln(k)-factln(n-k))); 

} 

double  factln(int  n) 

{  double  gammln(double  xx); 

static  double  a[101]; 
if  (n  <=  1)  return  0.0; 

if  (n  <=  100)  return  a[n]  ?  a[n] ;  (a[n]=gaininln(n+1.0)); 
else  return  ganmiln(n+1.0); 

} 

double  gammlnfdouble  xx) 

{  double  x,y,tmp,ser; 

static  double  cof[6]={76.18009172947146,-86.50532032941677, 
24.01409824083091,-1.231739572450155, 
0.1208650973866179e-2,-0.5395239384953e-5}; 

intj; 

y=x=xx; 

tnip=x+5.5; 

trap  -=  (x+0.5)*log(tmp); 

ser=1.000000000190015; 

for  (j=0y<=5  j++)  ser  +=  coflj]/++y; 

return  -tnip+^log(2.50662827463 10005*ser/x); 


#undefEPS 

#undefFPMIN 

#undefMAXIT 

#undefXMIN 

#undefPI 

#undefNRANSI 

#undefNUSEl 

#undefNUSE2 

#undefRTPI02 
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